

      program main

      use globais
      implicit none

!	inversao 
      double precision, dimension(:), allocatable :: X !< variaveis de otimizacao/parametros da inversao

      double precision:: cmin0,cmax0
      integer :: m !< numero de medidas 
      integer :: i,j
      double precision, dimension(:), allocatable :: Xaux !< variaveis de otimizacao/parametros da inversao auxiliar.


      double precision:: dreceiv !< distância entre receptores, variável auxiliar na geração dos sismogramas
      double precision, dimension(:), allocatable :: Gaux
      double precision :: Faux

      double precision :: VALUE	
      integer :: kncl
      integer :: knclbest !< Melhor número de classes para clusterização.
      doubleprecision :: valuebest !< Melhor valor da função objetivo entre os números de classe testados



	integer :: MaxIter

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
! ! non linear conjugate gradient, napack, netlib
! 	!input
	double precision :: STEP  !--STEP SIZE ALONG SEARCH DIRECTION FOR FINAL ITERATION		   
 	double precision :: Tol  !--COMPUTING TOLERANCE (ITERATIONS STOP WHEN MAX-NORM OF GRADIENT .LE. Tol)  
	integer :: limit
	double precision , dimension(:,:), allocatable::Work
	!output
	double precision :: E  !--MAX-NORM OF GRADIENT

! 
! 	EXTERNAL BOTH,GRAD,PRE,VALUE
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 


      !< L-BFGS 
!       PARAMETER(NDIM=2000,MSAVE=7,NWORK=NDIM*(2*MSAVE +1)+2*MSAVE)
      DOUBLE PRECISION, dimension(:), allocatable:: DIAG
      DOUBLE PRECISION, dimension(:), allocatable:: W
      DOUBLE PRECISION :: EPS,XTOL,GTOL,STPMIN,STPMAX
      INTEGER ::  IPRINT(2),IFLAGlbfgs,ICALL,MP,LP
      LOGICAL ::  DIAGCO
      integer :: msave
      integer :: icallmax


      EXTERNAL LB2
      COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 








      double precision :: lbeg,lend,ldnum

      double precision :: t ! provisorio     
       integer :: it ! provisorio 
       double precision :: dt ! para sismograma e snapshots     
      integer :: nt ! numero de instantes para sismograma 
      integer :: simulacao
      character*5 filename1,filename2,filename3
      character*1 :: barra= achar(92),pl= achar(39),plpl= achar(34)     
      character*2 :: charit

      double precision :: mapCemX

      complex*16 ::Ptrecaux
      double precision :: residuo
      character*5 :: chks
      character*2 :: chigsr


      call entrada(simulacao)


!       write(*,*)'# transformada de Fourier da fonte'
      call fonte()



	  if(simulacao.eq.0)then

	      write(*,*)'Solver Inversao'


	      call entradainversao(maxIter,step,tol,limit,cmin0,cmax0)


	      write(*,*)'alocando X...'
	      allocate(x(n))
	      write(*,*)'alocando Xaux...'
	      allocate(xaux(n))
	      write(*,*)'alocando cX...'
	      allocate(CX(n))

	      do j=1,nz0
		  do i=1,nx0
		      X((j-1)*nx0+i) = cmin0  + (cmax0-cmin0)/dfloat(nz0-1)*(j-1)
		  enddo
	      enddo
	      do i=1,N
		  X(i)    =  (     X(i) - (cmax+cmin)/2  )*2.2/(cmax-cmin)!isso eh uma aproximacao da funcao nao linear de X em C por uma reta
	      enddo


	      filename3= 'sol_c'
	      call system('rm '//filename3//'.bin '//filename3//'.ps')
	      open(unitsol,file=filename3//'.bin',status="unknown",form="unformatted",access="direct",recl=4*Nx0*Nz0)      

	      open(unitF,file='outfuncao.txt',status="unknown")    
	      open(unitFgrad,file='outfuncao_grad.txt',status="unknown") 


	      write(*,*)'alocando Xn...'
	      allocate(Xn(n))
	      write(*,*)'alocando cxn...'
	      allocate(CXn(n))
	      Xn = X



	      do kloop=1,Nloops

!  	      do kw=kloop,Nomega,Nloops !loop seletivo, seleciona freqs de acordo com o numero de loops
 	      do kw=1,Nomega
! 		    write(*,*) '#menor comprimento de onda (pelo cmin)= ', cmin/(omegas(kw)/(2*Pi))
! 		    write(*,*) '#intervalos por comprimento de onda (pelo cmin)= ',  ( cmin/(omegas(kw)/(2*Pi)))/delta
		    write(*,*)' kw ',kw,': ',omegas(kw)/(2.d0*Pi),'Hz'

		    norm2PBarra=0.d0
		    do igsr=1, gruposSR
			do ks=1,Nsource(igsr)
			    do kr=1,Nrec(igsr)
				norm2PBarra= norm2PBarra + PreceivBARRA(kr,ks,kw,igsr)*CONJG(PreceivBARRA(kr,ks,kw,igsr))
			    enddo
			enddo
		    enddo
		    write(*,*)'||Pbarra||**2=',norm2PBarra,'denominador de Fobj'


! 		    imetodo = 1
! 		    write(*,*)'alocando work para o cg do napack...'
! 		    allocate(Work(N,3))	
! 		    iterant = -1! numero da iteracao anterior
! 		    call CG(X,E,ITER,STEP,Tol,LIMIT,N,MaxIter,VALUE,GRAD,BOTH,PRE,Work)
! 		    write(*,*)'fim do CG   step:',step
! ! 		    step=0.0d0
! 		    deallocate(Work)! do cg do napack


		    imetodo = 2
		    iterant = -1! numero da iteracao anterior
		    MSAVE = 7
		    allocate(DIAG(N))
		    allocate(W(N*(2*MSAVE +1)+2*MSAVE))
		    IPRINT(1)= 1
		    IPRINT(2)= 0

		    DIAGCO= .FALSE.
		    EPS= 1.0D-5
		    XTOL= tol
		    ICALL=0
		    icallmax=25
		    maxiterLBFGS = limit
		    IFLAGlbfgs=0
		    DO
 			if(iflaglbfgs.eq.1 .or. (iflaglbfgs.eq.0 .and. diagco .eqv. .false.))then
			    call BOTH(fobj,gradFobj,X)
			else if(iflaglbfgs.eq.2 .or. (iflaglbfgs.eq.0 .and. diagco .eqv. .true.))then
			    call BOTHandH(fobj,gradFobj,X,DIAG)
			endif


			call LBFGS(N,Msave,iter,X,Fobj,GradFobj,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAGlbfgs)

			IF(IFLAGlbfgs.LE.0) exit

			ICALL=ICALL + 1

			write(*,*) 'maxiter',maxiterLBFGS,'iter',iter,'icall',icall
! 			IF(ICALL.GT.icallmax) exit ! We allow at most 2000 evaluations of F and G
			IF(iter.GT.maxiterLBFGS) exit ! We allow at most maxiterLBFGS iterations
		    ENDDO
		    deallocate(DIAG)
		    deallocate(W)


		    call BOTH(fobj,gradFobj,X)!<só para plotar a solucao final

		    open(unitX,file='outX.txt',status="unknown")
		    write(unitX,*)'# kw: ',kw,': ',omegas(kw)/(2.d0*Pi),'Hz',' X   final: '
		    do i=1,nx0*nz0
			write(unitX,*)X(i)
		    enddo
		    close(unitX)

		    call systempsmovieSOL(nx,nz,delta,delta,0.0d0,-dfloat(Npml-1)*delta,0.0d0,1.d0,&
			& filename3,'                    ','profundidade (m)    ',&
			&'distancia (m)       ',wbox,Lz0/4.d0,Lx0/4.d0,0.0d0,Lz0,0.0d0,Lx0,cmin,cmax) !figura c do pd sem PML

		    call systempsmovieSOL(nx,nz,delta,delta,0.0d0,-dfloat(Npml-1)*delta,0.0d0,1.d0,&
			& filename3,'                    ','profundidade (m)    ',&
			&'distancia (m)       ',wbox,Lz0/4.d0,Lx0/4.d0,0.0d0,&
			& Lz0+dfloat(Npml-1)*delta,-dfloat(Npml-1)*delta,Lx0+dfloat(Npml-1)*delta,cmin,cmax) !figura c do PD com PML

		    call systempsmovieSOL(nx0,nz0,deltax0,deltaz0,deltaz0/2.0d0,deltax0/2.0d0,0.0d0,1.d0,&
			& filename3,'                    ','profundidade (m)    ',&
			&'distancia (m)       ',wbox,Lz0/4.d0,Lx0/4.d0,0.0d0,Lz0,0.0d0,Lx0,cmin,cmax) !figura c0 do pi O(sem PML)




	      enddo ! loop das frequencias

			kw = kw-1

			call clusterizainterat(c0,nx0,nz0)
			call suaviza(c0,nx0,nz0)


! 			write(4,*)! salta linha no arquivo fort.4
! 			! open(1,file='velX0.bin',status="unknown",form="unformatted",access="direct",recl=4*Nx0*Nz0)
! 			! read(1,rec=1)((creal4(i,j),j=1,Nz0),i=1,Nx0)
! 			! close(1)
! 			! do j=1,nz0
! 			!     do i=1,nx0
! 			! 	c0(i,j) = creal4(i,j)
! 			!     enddo
! 			! enddo
! 			! do j=1,nz0
! 			!     do i=1,nx0
! 			! 	X((j-1)*nx0+i) = mapCemX(c0(i,j),cmin,cmax)
! 			!     enddo
! 			! enddo
! 			write(4,*)0,VALUE(X)
! 			knclbest = 0
! 			valuebest = 1000.d0
! 
! 			do kncl = 2,int(nz0/3) !< busca o melhor numero de classes para a clusterização.
! 			    call modelo2d(X)
! 			    ! open(1,file='velX0.bin',status="unknown",form="unformatted",access="direct",recl=4*Nx0*Nz0)
! 			    ! read(1,rec=1)((creal4(i,j),j=1,Nz0),i=1,Nx0)
! 			    ! close(1)
! 			    ! do j=1,nz0
! 			    !     do i=1,nx0
! 			    ! 	c0(i,j) = creal4(i,j)
! 			    !     enddo
! 			    ! enddo
! 			    call clusteriza(c0,nx0,nz0,kncl)
! ! 			    call suaviza(c0,nx0,nz0)
! 			    do j=1,nz0
! 				do i=1,nx0
! 				    Xaux((j-1)*nx0+i) = mapCemX(c0(i,j),cmin,cmax)
! 				enddo
! 			    enddo
! 			    write(4,*)'****************************',kncl,VALUE(Xaux)
! 			    if(VALUE(Xaux).lt.valuebest)then
! 				knclbest = kncl
! 				valuebest = VALUE(Xaux)
! 			    endif
! 			enddo
! 
! 			if(knclbest.gt.0)then !< se knclbest for igual a zero, não se faz clusterização e X não muda
! 			    call modelo2d(X)
! 			    call clusteriza(c0,nx0,nz0,knclbest)
! ! 			    call suaviza(c0,nx0,nz0)
! 			    do j=1,nz0
! 				do i=1,nx0
! 				    X((j-1)*nx0+i) = mapCemX(c0(i,j),cmin,cmax)
! 				enddo
! 			    enddo
! 			endif
! 			call BOTH(fobj,gradFobj,X)!<só para plotar a solucao final
! 			read(*,*)

	      enddo ! loop do numero de loops



! !		    PLOTAR FIGURA FINAL
! 		    iterant = -1! so p garantir a plotagem com esta reavaliacao final
! 		    write(*,*)' kw: ',kw,' X: ',X ! ta errado este kw, sai do loop valendo +1
! 		    allocate(Gaux(nx0*nz0))
! 		    call BOTH(Faux,Gaux,X) !so para salvar a fobj do resultado final
! 		    deallocate(Gaux)
! 		    call systempsmovieSOL(nx,nz,delta,0.0d0,-dfloat(Npml-1)*delta,0.0d0,1.d0,&
! 			& 'sol_c','                    ','profundidade (m)    ',&
! 			&'distancia (m)       ',wbox,Lz0/4.d0,Lx0/4.d0,cmin,cmax)


! 		    call refinagrid(X,n,m) ! entra n antigo e volta recalculado
! 		    deallocate(X)
! 		    n = nx0* nz0
! 		    allocate(X(n))  ! aloca novo X e reinicia as variaveis
! 		    do j=1,nz0
! 			do i=1,nx0
! 			    X((j-1)*nx0 +i) = c0(i,j)
! 			enddo
! 		    enddo
! 		    ftol = ftol/2.
! 		    xtol = xtol/2.
! 		    gtol = gtol/2.
! 		    maxfev = 10*n+n



	      deallocate(c0)


	      deallocate(X)
	      deallocate(Xaux)
	      deallocate(CX)
	      deallocate(xref)
	      deallocate(CXR)
	      deallocate(Xn)
	      deallocate(CXn)
	      deallocate(resid)
! 	      deallocate(Jacob)
! 	      deallocate(Hprecond)

	      deallocate(deltaR) ! matriz: cada coluna é vetor delta correspondente a um receptor      
	      deallocate(JdeX)     ! du/dX
	      deallocate(residexpandido)
	      deallocate(GradFobj) 
	      deallocate(GradTreg)
	      deallocate(depX)
	      deallocate(PreceivBARRA) 
	      deallocate(Hvalue)
	      deallocate(Hrowind)
	      deallocate(Hcolind)

	      deallocate(f) !lado direito do sistema do problema direto
	      deallocate(rowind)!CCS da matriz de impedancia
	      deallocate(colptr)!CCS da matriz de impedancia
	      deallocate(values)!valores nao nulos da matriz do sistema do problema direto

   
	      close(unitsol)
	      close(unitF)  
	      close(unitFgrad)




	  elseif(simulacao.eq.1)then


	      write(*,*) '# simulacao'

! 	      call modelo2dhustedt()
	      call modelo2dcamadas()

	      write(*,*) '#menor comprimento de onda (no modelo) = ', cminreal/(omegas(Nomega)/(2*Pi))
	      write(*,*) '#intervalos por comprimento de onda (no modelo)= ',  ( cminreal/(omegas(Nomega)/(2*Pi)))/delta

	      call systempsimage(c,nx,nz,delta,delta,0.d0,-(Npml-1)*delta,'veloc',&
		& 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
		& wbox,Lz0/4.d0,Lx0/4.d0, 0.0d0,Lz0,0.0d0,Lx0)!modelo do pd sem PML

	      call systempsimage(rho,nx,nz,delta,delta,0.d0,-(Npml-1)*delta,'densi',&
		& 'velocidade (m/s)    ','profundidade (m)    ', 'distancia (m)       ',&
		& wbox,Lz0/4.d0,Lx0/4.d0, 0.0d0,Lz0,0.0d0,Lx0)!modelo do pd sem PML



! 	      write(*,*)'alocando P(Nx,Nz,Nomega,NsourceSnap)'
! 	      allocate(P(Nx,Nz,Nomega,NsourceSnap))

	      write(*,*)'alocando P(Nx,Nz,NsourceSnap,gruposSR)'
	      allocate(P(Nx,Nz,NsourceSnap,gruposSR))
	      write(*,*)'alocando Ptreal(Nx,Nz,Ntsnaps,NsourceSnap,gruposSR)'
	      allocate(Ptreal(Nx,Nz,Ntsnaps,NsourceSnap,gruposSR))
! 	      allocate(Ptimag(Nx,Nz,Ntsnaps,NsourceSnap))
! 	      dtsnaps = (tf/dfloat(ntsnaps-1))

! 	      ideltasnap = int(Nsource/(Nsourcesnap-1))


	      nt = 2*Nomega+1
	      dt=(tf/dfloat(nt))

	      do kw = 1,Nomega

		    write(*,*)'.... kw',kw
		    call pd2d(omegas(kw),0,0,0)!< iflagF=0 iflagG=0 iflagH=0

		    do igsr=1,gruposSR !loop nos grupos fonte receptor
			ideltasnap = int(Nsource(igsr)/(Nsourcesnap-1))
			do kss=1,NsourceSnap ! acumula resp no dom do tempo p cada freq resolvida para algumas fontes
			    ks=(kss-1)*ideltaSnap + 1
			    if(ks.gt.Nsource(igsr))ks=Nsource(igsr)
			    do it=1,ntsnaps
    ! 			    t = (it-1)*dtsnaps
				t=(it-1)*nt/(ntsnaps-1)*dt ! para alguns ntsnaps instantes usados na transformada da fonte 
				do j=1,Nz
				    do i=1,Nx
					! parte real do campo depois da tranformada inversa
					Ptreal(i,j,it,kss,igsr) = Ptreal(i,j,it,kss,igsr) + real(          P(i,j,kss,igsr)*exp(imag*omegas(kw)*t) + &  
											& CONJG(P(i,j,kss,igsr))*exp(imag*omegas(2*Nomega-(kw-1))*t))
					! parte imaginaria do campo (praticamente nula ja que f eh real)
    ! 				    Ptimag(i,j,it,kss) = Ptimag(i,j,it,kss) +aimag(  P(i,j,kss)*exp(imag*omegas(kw)*t) + &  
    ! 										    & CONJG(P(i,j,kss))*exp(imag*omegas(2*Nomega-(kw-1))*t)) 
				    enddo! em x
				enddo! em z
			    enddo ! fim do loop nos instantes
			enddo! nas fontes escolhidas para os snapshots
		    enddo ! loop os grupos fonte receptor


	      enddo! loop das frequencias


	      open(1,file='medidasfreq.txt')
	      write(1,*)'#Nreceptor Nfonte Nfrequencia Medida'
	      do kw=1,Nomega
		  do igsr=1,gruposSR
		      do ks=1,Nsource(igsr)
			  do kr=1,Nrec(igsr)
			      write(1,*)Preceiv(kr,ks,kw,igsr)
    ! 			  write(1,101)kr,ks,kw,Preceiv(kr,ks,kw)
			  enddo
		      enddo
		  enddo
	      enddo
	      close(1)

	      ! Faz a transformada inversa calculando a resp no dominio do tempo em cada receptor e gera sismograma
	      write(*,*)'FFT inversa e sismogramas'

	      nt = 2*Nomega+1
	      dt=(tf/dfloat(nt))


	      do igsr=1,gruposSR

		  allocate(PtrecR(nt,Nrec(igsr)))
		  allocate(PtrecI(nt,Nrec(igsr)))
		  ideltasnap = int(Nsource(igsr)/(Nsourcesnap-1))

		  do kss=1,NsourceSnap

		      ks=(kss-1)*ideltaSnap + 1
		      if(ks.gt.Nsource(igsr))ks=Nsource(igsr)

		      do it=1,nt
			  t = (it-1)*dt
			  do kr=1,Nrec(igsr)
			      Ptrecaux = 0.0
			      do kw = 1,Nomega
				      Ptrecaux = Ptrecaux + Preceiv(kr,ks,kw,igsr)*exp(imag*omegas(kw)*t)+&
					    & conjg(Preceiv(kr,ks,kw,igsr))*exp(imag*omegas(2*Nomega-(kw-1))*t)
			      enddo
			      PtrecR(it,kr) = real(Ptrecaux)
			      PtrecI(it,kr) =aimag(Ptrecaux)

			  enddo
		      enddo

		      if(ks.lt.10)then
			  write(chks,'(i1)')ks
			  chks = '0000'//chks
		      elseif(ks.lt.100)then
			  write(chks,'(i2)')ks
			  chks = '000'//chks
		      elseif(ks.lt.1000)then
			  write(chks,'(i3)')ks
			  chks = '00'//chks
		      elseif(ks.lt.10000)then
			  write(chigsr,'(i4)')ks
			  chks = '0'//chks
		      elseif(ks.lt.100000)then
			  write(chks,'(i5)')ks
		      endif
		      if(igsr.lt.10)then
			  write(chigsr,'(i1)')igsr
			  chigsr = '0'//chigsr
		      elseif(igsr.lt.100)then
			  write(chigsr,'(i2)')igsr
		      endif

		      dreceiv=dfloat(ireceiv(2,igsr)-ireceiv(1,igsr))*delta ! distancia entre receptores
		      call systempswigp(PtrecR,nt,nrec,dt,dreceiv,&
				& 0.0d0,dfloat(ireceiv(1,igsr)-isource(ks,igsr))*delta,&
				& 'sismR'//chigsr//chks,'      sismograma     ','      tempo (s)      ','     offset (m)      ',&
				& wbox,tf/5.0d0,dfloat(ireceiv(nrec(igsr),igsr)-ireceiv(1,igsr))*delta/5.d0)


		  enddo! loop das fontes nsourcesnap
		  deallocate(PtrecR)
		  deallocate(PtrecI)

	      enddo ! fim do loop dos grupos fonte receptor
	      !  fim da geração dos sismogramas


	      do igsr=1,gruposSR
	      ! a partir da resp ja no dom do tempo  gera os snapshots
		  write(*,*)'snapshots'
		  ideltasnap = int(Nsource(igsr)/(Nsourcesnap-1))
		  if(igsr.lt.10)then
		      write(chigsr,'(i1)')igsr
		      chigsr = '0'//chigsr
		  elseif(igsr.lt.100)then
		      write(chigsr,'(i2)')igsr
		  endif
		  call systempsmovie(Ptreal(1,1,1,1,igsr),nx,nz,nomega,Nsource,delta,0.00d0,-dfloat(Npml-1)*delta,0.0d0,tf,&
					  &'        campo       ','  profundidade (m)  ',&
					  &'   distancia (m)    ',wbox,Lz0/4.d0, Lx0/4.d0 ,  &
					  & ideltasnap, nsourcesnap, ntsnaps,'snapR'//chigsr  )
		  !   fim da geracao dos snapshots para algumas fontes
	      enddo

	      deallocate(P)
	      deallocate(Ptreal)
! 	      deallocate(Ptimag)


	  endif !fim do if da simulacao	  


      deallocate(csix)
      deallocate(csiz)
      deallocate(c)

      deallocate(creal4)

      deallocate(rho)
      deallocate(b)
      deallocate(K)





      deallocate(Preceiv)





      deallocate(Sw)
      deallocate(omegas)
      deallocate(ireceiv)
      deallocate(jreceiv)






      stop

 101  format(3i10,'  (',E15.6,' , ',E15.6,' ) ')
      end program main